/* TridiagonalSystem.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#ifndef _TridiagonalSolver_
#define _TridiagonalSolver_
#include <src/miscellaneous/OS_basic.h>

/*! \brief Tri Diagonal solver iterator oriented

	It is a iterator templated version of a classical tridiagonal solver using iterator definition to solve the tridiagonal system
	Au=b.
	
*/

template <class data, class Container>
void TridiagonalSolverIterator(
							   typename data::iterator 				result_iterator, 
							   typename data::const_iterator 		rhs_iterator, 
							   typename Container::const_iterator 	diagonal_iterator, 
							   typename Container::const_iterator 	upper_diagonal_iterator, 
							   typename Container::const_iterator 	lower_diagonal_iterator, 
							   OS_size						size)
{
		typedef Container::const_iterator	const_iterator;
		typedef Container::iterator			iterator;
		typedef Container::value_type		value_type;
		typedef data::const_iterator		data_const_iterator;
		typedef data::iterator				data_iterator;

        //take care : no run time check here

		Container tmp(size);

		const_iterator iterator_end(diagonal_iterator);
		iterator_end += size;
		iterator tmp_iterator = tmp.begin();

		data_iterator result_iterator_1(result_iterator);

        value_type bet		= (*diagonal_iterator);
		(*result_iterator)	= (*rhs_iterator)/bet;

        for (++diagonal_iterator
			, ++tmp_iterator
			, ++result_iterator
			,++rhs_iterator
			; diagonal_iterator!=iterator_end
			; ++diagonal_iterator
			, ++tmp_iterator
			, ++result_iterator
			, ++rhs_iterator
			, ++upper_diagonal_iterator
			, ++lower_diagonal_iterator
			, ++result_iterator_1)
		{
            (*tmp_iterator)		= (*upper_diagonal_iterator)/bet;
            bet					= (*diagonal_iterator)-(*lower_diagonal_iterator)*(*tmp_iterator);
            (*result_iterator)	= ( (*rhs_iterator)-(*lower_diagonal_iterator)*(*result_iterator_1) )/bet;
        }
		--result_iterator;
		result_iterator_1 = result_iterator;
		iterator_end = tmp_iterator;
		iterator_end -= size;

        for (--result_iterator,
			-- tmp_iterator
			; tmp_iterator != iterator_end
			; --tmp_iterator
			, --result_iterator
			, --result_iterator_1) 
		{
            (*result_iterator) -= (*tmp_iterator)*(*result_iterator_1);
		}
};



/*! \brief Tri Diagonal solver

	It is a iterator templated version of a classical tridiagonal solver using iterator definition to solve the tridiagonal system
	Au=b. It uses the proxy and call TridiagonalSolverIterator.
	
*/
template <class data, class Container>
data TridiagonalSolver(const data & rhs, const Container & diagonal, const Container & upper_diagonal,const Container & lower_diagonal)
{

		data result(rhs.size_list());
		OS_size size = diagonal.size();


		TridiagonalSolverIterator<data,Container>(result.begin(),rhs.begin(),diagonal.begin(),upper_diagonal.begin(),lower_diagonal.begin(),size);

		return result;

};
/*! \brief Tri Diagonal matrix action

	It is a templated algorithm using iterator definition to compute
	Au=b, where A is a tridiagonal matrix (upper_diagonal, diagonal, lower_diagonal) and u is a vector (rhs).
	No run time check performed, thus take care of the input.
*/
template <class data, class Container>
data Apply(const data & rhs, const Container & diagonal, const Container & upper_diagonal,const Container & lower_diagonal)
{

        //take care : no run time check here
		data result(rhs.size_list());
		OS_size size = diagonal.size();

		typedef Container::const_iterator	const_iterator;
		typedef Container::iterator			iterator;
		typedef Container::value_type		value_type;
		typedef data::const_iterator		data_const_iterator;
		typedef data::iterator				data_iterator;

		const_iterator iterator_				= diagonal.begin(); 
		const_iterator iterator_end				= diagonal.end();

		data_const_iterator rhs_iterator = rhs.begin();
		data_iterator result_iterator = result.begin();


        for (
			; iterator_!=iterator_end
			; ++iterator_
			, ++rhs_iterator
			, ++result_iterator)
		{
            *result_iterator = (*iterator_)*(*rhs_iterator);
        }

		iterator_							= upper_diagonal.begin(); 
		iterator_end						= iterator_, iterator_end += (size-1);
		rhs_iterator						= rhs.begin(), ++rhs_iterator;
		result_iterator						= result.begin();
        for (
			; iterator_!=iterator_end
			; ++iterator_
			, ++rhs_iterator
			, ++result_iterator)
		{
            *result_iterator += (*iterator_)*(*rhs_iterator);
        }

		iterator_							= lower_diagonal.begin(); 
		iterator_end						= iterator_, iterator_end += (size-1);
		rhs_iterator						= rhs.begin();
		result_iterator						= result.begin(), ++result_iterator;
        for (
			; iterator_!=iterator_end
			; ++iterator_
			, ++rhs_iterator
			, ++result_iterator)
		{
            *result_iterator += (*iterator_)*(*rhs_iterator);
        }

        return result;
};


#endif